Load Data

Single-cell data

set.seed(1234)

# Load Data
WT_Mouse <- readRDS("./01.2_main_E17andPN_WT_AddMeta_Annotated.RDS")

WT_Mouse$NewCelltype <- ""
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 0)] <- "(0) E ADR Chromaffin"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 1)] <- "(1) PN NOR Chromaffin"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 2)] <- "(2) PN ADR Chromaffin"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 3)] <- "(3) Endothelial"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 4)] <- "(4) Macrophage"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 5)] <- "(5) E NOR Chromaffin"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 6)] <- "(6) E Glia|SCP"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 7)] <- "(7) Cortex ZG"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 8)] <- "(8) Cycling TH+"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 9)] <- "(9) Cortex X Zone"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 10)] <- "(10) E Neuroblast"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 11)] <- "(11) PN Sustentacular"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 12)] <- "(12) Cortex ZF"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 13)] <- "(13) Macrophage"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 14)] <- "(14) Cortex ZG"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 15)] <- "(15) Dendritic cell"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 16)] <- "(16) T cell"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 17)] <- "(17) NK cell"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 18)] <- "(18) Mesenchymal"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 19)] <- "(19) Dendritic cell"

WT_Mouse$DetailedAge <- ""
WT_Mouse$DetailedAge[which(WT_Mouse$wenyu_names %in%c ("E17nTAM_11_5", "E17nTAM_8_5"))] <- "Embryonic (E17)"
WT_Mouse$DetailedAge[which(WT_Mouse$wenyu_names %in%c ("P5"))] <- "Postnatal (P5)"
WT_Mouse$DetailedAge[which(WT_Mouse$wenyu_names %in%c ("P14"))] <- "Postnatal (P14)"
WT_Mouse$DetailedAge[which(WT_Mouse$wenyu_names %in%c ("WT_AG_3m"))] <- "Postnatal (P90)"
WT_Mouse$DetailedAge[which(WT_Mouse$wenyu_names %in%c ("AdultAG"))] <- "Postnatal (Aged)"
WT_Mouse$DetailedAge <- factor(WT_Mouse$DetailedAge, levels = c("Embryonic (E17)", "Postnatal (P5)", "Postnatal (P14)", "Postnatal (P90)", "Postnatal (Aged)"))

WT_Mouse$GrossAge <- ""
WT_Mouse$GrossAge[which(WT_Mouse$wenyu_names %in%c ("E17nTAM_11_5", "E17nTAM_8_5"))] <- "Embryonic (E17)"
WT_Mouse$GrossAge[which(WT_Mouse$wenyu_names %in%c ("P5", "P14", "WT_AG_3m", "AdultAG"))] <- "Postnatal"
WT_Mouse$GrossAge <- factor(WT_Mouse$GrossAge, levels = c("Embryonic (E17)", "Postnatal"))

NewCelltypes <- sort(unique(WT_Mouse$NewCelltype))
NewCelltypesOrdered <- NewCelltypes[c(1,2,13:20,3:12)]

WT_Mouse$NewCelltype <- factor(WT_Mouse$NewCelltype, levels = NewCelltypesOrdered)
WT_Mouse$seurat_clusters <- factor(WT_Mouse$seurat_clusters, levels = c(0:19))

Idents(WT_Mouse) <- 'NewCelltype'
DefaultAssay(WT_Mouse) <- 'RNA'
WT_Mouse
## An object of class Seurat 
## 38268 features across 7101 samples within 1 assay 
## Active assay: RNA (38268 features, 2500 variable features)
##  3 layers present: counts, data, scale.data
##  2 dimensional reductions calculated: pca, umap
a <- DimPlot(WT_Mouse, reduction = "umap", group.by = "NewCelltype")
pbuild <- ggplot2::ggplot_build(a)
pdata <- pbuild$data[[1]]
pdata <-  pdata[order(pdata$group), ]
ucols <- unique(pdata$colour)
names(ucols) <- levels(WT_Mouse$NewCelltype)
ucols
##  (0) E ADR Chromaffin (1) PN NOR Chromaffin (2) PN ADR Chromaffin 
##             "#F8766D"             "#EA8331"             "#D89000" 
##       (3) Endothelial        (4) Macrophage  (5) E NOR Chromaffin 
##             "#C09B00"             "#A3A500"             "#7CAE00" 
##        (6) E Glia|SCP         (7) Cortex ZG       (8) Cycling TH+ 
##             "#39B600"             "#00BB4E"             "#00BF7D" 
##     (9) Cortex X Zone     (10) E Neuroblast (11) PN Sustentacular 
##             "#00C1A3"             "#00BFC4"             "#00BAE0" 
##        (12) Cortex ZF       (13) Macrophage        (14) Cortex ZG 
##             "#00B0F6"             "#35A2FF"             "#9590FF" 
##   (15) Dendritic cell           (16) T cell          (17) NK cell 
##             "#C77CFF"             "#E76BF3"             "#FA62DB" 
##      (18) Mesenchymal   (19) Dendritic cell 
##             "#FF62BC"             "#FF6A98"

Cell-cell-communication All Cell types

CellChat_AgedAdult <- readRDS("./CellChat_CellCellCommunication/WT_Mouse_AllCellType_SpeAge/WT_Mouse_AgedAdult_AllCelltype_CCC.rds")
CellChat_P90 <- readRDS("./CellChat_CellCellCommunication/WT_Mouse_AllCellType_SpeAge/WT_Mouse_3MonthAdult_AllCelltype_CCC.rds")
CellChat_P14 <- readRDS("./CellChat_CellCellCommunication/WT_Mouse_AllCellType_SpeAge/WT_Mouse_PN_14D_AllCelltype_CCC.rds")
CellChat_P5 <- readRDS("./CellChat_CellCellCommunication/WT_Mouse_AllCellType_SpeAge/WT_Mouse_PN_5D_AllCelltype_CCC.rds")
CellChat_E17 <- readRDS("./CellChat_CellCellCommunication/WT_Mouse_AllCellType_SpeAge/WT_Mouse_E_17_AllCelltype_CCC.rds")
CellChat_AgedAdult@meta$NewCelltype <- as.character(WT_Mouse@meta.data[rownames(CellChat_AgedAdult@meta), "NewCelltype"])
CellChat_AgedAdult <- setIdent(CellChat_AgedAdult, ident.use = "NewCelltype")
## Rename cell groups but do not change the order!
## The cell group order before renaming is  1-PN_NOR Chromaffin 11-PN_Glia/Sustentacular 12-PN_Cortex ZF 15-PN_Dendritic cell 16-PN_T cell 17-PN_NK cell 18-PN_Mesenchymal 19-PN_Dendritic cell 2-PN_ADR Chromaffin 3-PN|E_Endothelial 4-PN_Macrophage 7-PN_Cortex ZG 
## The cell group order after renaming is  (1) PN NOR Chromaffin (11) PN Sustentacular (12) Cortex ZF (15) Dendritic cell (16) T cell (17) NK cell (18) Mesenchymal (19) Dendritic cell (2) PN ADR Chromaffin (3) Endothelial (4) Macrophage (7) Cortex ZG
## Warning in setIdent(CellChat_AgedAdult, ident.use = "NewCelltype"): All the calculations after `computeCommunProb` should be re-run!!
##     These include but not limited to `computeCommunProbPathway`,`aggregateNet`, and `netAnalysis_computeCentrality`.
CellChat_P90@meta$NewCelltype <- as.character(WT_Mouse@meta.data[rownames(CellChat_P90@meta), "NewCelltype"])
CellChat_P90 <- setIdent(CellChat_P90, ident.use = "NewCelltype")
## Rename cell groups but do not change the order!
## The cell group order before renaming is  1-PN_NOR Chromaffin 11-PN_Glia/Sustentacular 15-PN_Dendritic cell 16-PN_T cell 17-PN_NK cell 2-PN_ADR Chromaffin 3-PN|E_Endothelial 4-PN_Macrophage 
## The cell group order after renaming is  (1) PN NOR Chromaffin (11) PN Sustentacular (15) Dendritic cell (16) T cell (17) NK cell (2) PN ADR Chromaffin (3) Endothelial (4) Macrophage
## Warning in setIdent(CellChat_P90, ident.use = "NewCelltype"): All the calculations after `computeCommunProb` should be re-run!!
##     These include but not limited to `computeCommunProbPathway`,`aggregateNet`, and `netAnalysis_computeCentrality`.
CellChat_P14@meta$NewCelltype <- as.character(WT_Mouse@meta.data[rownames(CellChat_P14@meta), "NewCelltype"])
CellChat_P14 <- setIdent(CellChat_P14, ident.use = "NewCelltype")
## Rename cell groups but do not change the order!
## The cell group order before renaming is  1-PN_NOR Chromaffin 11-PN_Glia/Sustentacular 13-PN_Macrophage 14-PN_Cortex ZG 15-PN_Dendritic cell 16-PN_T cell 2-PN_ADR Chromaffin 3-PN|E_Endothelial 4-PN_Macrophage 7-PN_Cortex ZG 8-E_Cycling Medullary cell 9-PN_Cortex X Zone 
## The cell group order after renaming is  (1) PN NOR Chromaffin (11) PN Sustentacular (13) Macrophage (14) Cortex ZG (15) Dendritic cell (16) T cell (2) PN ADR Chromaffin (3) Endothelial (4) Macrophage (7) Cortex ZG (8) Cycling TH+ (9) Cortex X Zone
## Warning in setIdent(CellChat_P14, ident.use = "NewCelltype"): All the calculations after `computeCommunProb` should be re-run!!
##     These include but not limited to `computeCommunProbPathway`,`aggregateNet`, and `netAnalysis_computeCentrality`.
CellChat_P5@meta$NewCelltype <- as.character(WT_Mouse@meta.data[rownames(CellChat_P5@meta), "NewCelltype"])
CellChat_P5 <- setIdent(CellChat_P5, ident.use = "NewCelltype")
## Rename cell groups but do not change the order!
## The cell group order before renaming is  1-PN_NOR Chromaffin 11-PN_Glia/Sustentacular 2-PN_ADR Chromaffin 8-E_Cycling Medullary cell 
## The cell group order after renaming is  (1) PN NOR Chromaffin (11) PN Sustentacular (2) PN ADR Chromaffin (8) Cycling TH+
## Warning in setIdent(CellChat_P5, ident.use = "NewCelltype"): All the calculations after `computeCommunProb` should be re-run!!
##     These include but not limited to `computeCommunProbPathway`,`aggregateNet`, and `netAnalysis_computeCentrality`.
CellChat_E17@meta$NewCelltype <- as.character(WT_Mouse@meta.data[rownames(CellChat_E17@meta), "NewCelltype"])
CellChat_E17 <- setIdent(CellChat_E17, ident.use = "NewCelltype")
## Rename cell groups but do not change the order!
## The cell group order before renaming is  0-E_ADR Chromaffin 10-E_Neuroblast 3-PN|E_Endothelial 5-E_NOR Chromaffin 6-E_Glia/SCP 7-PN_Cortex ZG 8-E_Cycling Medullary cell 
## The cell group order after renaming is  (0) E ADR Chromaffin (10) E Neuroblast (3) Endothelial (5) E NOR Chromaffin (6) E Glia|SCP (7) Cortex ZG (8) Cycling TH+
## Warning in setIdent(CellChat_E17, ident.use = "NewCelltype"): All the calculations after `computeCommunProb` should be re-run!!
##     These include but not limited to `computeCommunProbPathway`,`aggregateNet`, and `netAnalysis_computeCentrality`.
CellChat_AgedAdult_WNT <- CellChat_AgedAdult
CellChat_AgedAdult_WNT@LR$LRsig <- CellChat_AgedAdult@LR$LRsig[which(CellChat_AgedAdult@LR$LRsig$ligand == "Wnt6"),]

CellChat_P90_WNT <- CellChat_P90
CellChat_P90_WNT@LR$LRsig <- CellChat_P90@LR$LRsig[which(CellChat_P90@LR$LRsig$ligand == "Wnt6"),]

CellChat_P14_WNT <- CellChat_P14
CellChat_P14_WNT@LR$LRsig <- CellChat_P14@LR$LRsig[which(CellChat_P14@LR$LRsig$ligand == "Wnt6"),]

CellChat_P5_WNT <- CellChat_P5
CellChat_P5_WNT@LR$LRsig <- CellChat_P5@LR$LRsig[which(CellChat_P5@LR$LRsig$ligand == "Wnt6"),]

CellChat_E17_WNT <- CellChat_E17
CellChat_E17_WNT@LR$LRsig <- CellChat_E17@LR$LRsig[which(CellChat_E17@LR$LRsig$ligand == "Wnt6"),]

Cell-cell-communication Medullary Cell types

cellChat_Merged <- readRDS("./CellChat_CellCellCommunication/WT_Mouse_MedullaCell_SpeAge/WT_Mouse_Merged_AllAgeGroup_CCC.Rds")

Figure 5 A

DimPlot(WT_Mouse, group.by='NewCelltype', reduction = 'umap', label = T, repel = T) + theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), axis.line = element_blank())

Figure 5 B

Medullary Cell types

WT_Mouse_Medulla <- subset(WT_Mouse, seurat_clusters %in% c(0, 1, 2, 5, 6, 8, 10, 11))
Medulla_CellNames <- colnames(WT_Mouse_Medulla)[which((2 * WT_Mouse_Medulla@reductions$umap@cell.embeddings[,1] + WT_Mouse_Medulla@reductions$umap@cell.embeddings[,2]) > 0)]
WT_Mouse_Medulla <- subset(WT_Mouse_Medulla, cells = Medulla_CellNames)

DimPlot(WT_Mouse_Medulla, group.by='GrossAge', pt.size=0.4, reduction = 'umap', label = F, repel = T)

Figure 5 C

OutDir <- "./DEG_GO/"
WT_Mouse_SpeDEGs <- SpecificDEG(WT_Mouse, Cluster = "NewCelltype", min.logfc = 0.5, p.adj = 0.05, slot = "data", save = T, return_res = T, outdir = OutDir, Prefix = "Whole_Data_All_Cluster")
Glia_DEGs <- FindMarkers(WT_Mouse, ident.1 = "(11) PN Sustentacular", ident.2 = "(6) E Glia|SCP", logfc.threshold = 0.5, group.by = "NewCelltype", only.pos = F)
saveRDS(Glia_DEGs, "Glia_Compararsion_PN_vs_E17_DEGs.rds")

head(Glia_DEGs)
# Define the pattern strings
patterns <- c("^\\(11\\) PN Sustentacular", "^\\(6\\) E Glia\\|SCP")

# Function to split specified columns and filter based on the patterns
split_and_extract <- function(df, columns_to_split, patterns, columns_to_keep) {
  # Initialize a copy of the dataframe to store results
  df_extracted <- df[,columns_to_keep, drop = FALSE]
  
  # Only split the specified columns
  for (col in columns_to_split) {
    df_extracted[, col] <- sapply(df[[col]], function(col_value) {
      # Split each column element based on "| " pattern
      col_split <- strsplit(col_value, "\\| ")
      
      # Filter the values that match the patterns
      filtered_values <- sapply(col_split, function(x) {
        matched <- x[grepl(paste(patterns, collapse = "|"), x)]
        if (length(matched) > 0) {
          return(matched)
        } else {
          return(NA)  # Return NA if no match
        }
      })
      filtered_values <- as.numeric(gsub(paste(paste0(patterns, ": "), collapse = "|"),"",filtered_values))
      return(filtered_values) 
    })
  }
  
  return(df_extracted)
}

# Specify the columns you want to split
columns_to_split <- c("others.pct.2", "others.pval", "others.padj", "others.log2FC")  # You can change this to your column names

# Specify the columns you want to keep
columns_to_keep <- c("gene", "cluster", "pct.1")  # You can change this to your column names

Genes2Show <- c("Sox10", "Wwtr1", "Tnc", "Mpz", "Lmo4", "Mef2c", "Gpr17", "Mmp2", "Ngfr", "Sox5", "Ctnnb1", "Moxd1", "Hif3a", "Atp1a2", "Olfml2a", "Col28a1", "Aspa", "Notch1", "Heyl", "Sfrp5", "Sfrp1", "Wnt6", "Gfap", "Tgfb2", "Sema3c", "Adam11", "Fam19a5", "Sparc", "Hey2", "Rgmb")

Glia_SepDEGs <- WT_Mouse_SpeDEGs[which(WT_Mouse_SpeDEGs$cluster %in% c("(11) PN Sustentacular", "(6) E Glia|SCP")),c("gene", "cluster", "pct.1", "others.pct.2", "others.pval", "others.padj", "others.log2FC")]

# Apply the function to extract the desired columns
result <- split_and_extract(Glia_SepDEGs, columns_to_split, patterns)
colnames(result) <- gsub("others.", "", colnames(result))
result$cluster2 <- ""
result$cluster2[which(result$cluster == "(6) E Glia|SCP")] <- "(11) PN Sustentacular"
result$cluster2[which(result$cluster == "(11) PN Sustentacular")] <- "(6) E Glia|SCP"

colnames(result)[2] <- "cluster1"
result <- result[,c(1,2,8,3:7)]
# View the result
# Create the volcano plot
result2 <- result
result2$log2FC[which(result2$cluster1 == "(6) E Glia|SCP")] <- -result2$log2FC[which(result2$cluster1 == "(6) E Glia|SCP")]
head(result2, 20)
Glia_DEGs_M <- Glia_DEGs[which(Glia_DEGs$p_val_adj != 1),]
colnames(Glia_DEGs_M) <- c("pval", "log2FC", "pct.1", "pct.2", "padj")
Glia_DEGs_M$gene <- rownames(Glia_DEGs_M)
Glia_DEGs_M$cluster1 <- "(11) PN Sustentacular"
Glia_DEGs_M$cluster2 <- "(6) E Glia|SCP"
Glia_DEGs_M <- Glia_DEGs_M[,colnames(result2)]
Glia_DEGs_M <- Glia_DEGs_M[-which(Glia_DEGs_M$gene %in% result2$gene),]
Glia_DEGs_M$DEG_Type <- "Non-Specific"
head(Glia_DEGs_M)
result2$DEG_Type <- "(6) E Glia|SCP"
result2$DEG_Type[which(result2$cluster1 == "(11) PN Sustentacular")] <- "(11) PN Sustentacular"
Glia_DEGs_Combined <- rbind(Glia_DEGs_M, result2)
Glia_DEGs_Combined$DEG_Type <- factor(Glia_DEGs_Combined$DEG_Type, levels = c("(11) PN Sustentacular", "(6) E Glia|SCP", "Non-Specific"))

# Create the volcano plot
p <- ggplot(Glia_DEGs_Combined, aes(x = log2FC, y = -log10(padj), color = DEG_Type)) + geom_point(alpha = 0.8, aes(size = ifelse(grepl("Glia-Specific", DEG_Type), 0.5, 0.2)), shape = 16) + 
  scale_color_manual(values = c("#01BBE1", "#3DB54A", "grey88")) + theme_minimal() + labs(title = "Volcano Plot", x = "log2 Fold Change", y = "-log10 Adjusted p-value") + 
  theme(legend.position = c(0.9, 0.92), plot.title = element_text(hjust = 0.5), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(color = "black"), 
        axis.ticks = element_line(color = "black"), legend.text = element_text(face="bold",colour = "black", size = 10), panel.border = element_rect(color = "black", fill = NA)) +
  geom_vline(xintercept = c(-0.5, 0.5), linetype = "dashed", color = "black") + geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "black") +
  geom_text_repel(data = result2[result2$gene %in% Genes2Show,], mapping = aes(x = log2FC, y = -log10(padj), color = DEG_Type, label = gene, size = 3), 
                  max.overlaps = Inf, show.legend = FALSE, fontface = "bold.italic") + 
  theme(strip.text.x = element_text(angle=0, face="bold",size = 15), strip.text.y=element_text(angle=90,face="bold",size = 30), title = element_text(face="bold",colour = "black", size = 12)) +
  theme(axis.text.x = element_text(face = "bold",colour = "black", angle = 0, hjust = 1, size = 10), axis.text.y = element_text(face="bold",colour = "black", size = 10)) + 
  guides(size = "none", color = guide_legend(override.aes = list(size = 4)))
p

Figure 5 D

color.heatmap <- "viridis"
color.use <- tryCatch({
  RColorBrewer::brewer.pal(n = 5, name = color.heatmap)
}, error = function(e) {
  (scales::viridis_pal(option = color.heatmap, direction = -1))(5)
})

MedullaCellAnnoOrder <- c("(6) E Glia|SCP", "(11) PN Sustentacular", "(10) E Neuroblast", "(8) Cycling TH+", "(0) E ADR Chromaffin", "(5) E NOR Chromaffin", "(2) PN ADR Chromaffin", "(1) PN NOR Chromaffin")
SpeDEG_Dotplot <- MyDotPlots(WT_Mouse_Medulla, result2[result2$gene %in% Genes2Show,],GeneCol = 'gene',Assay = "RNA", CelltypeCol = 'cluster1', Return = T,Draw = F, Group = "NewCelltype", CellTypeOrder = MedullaCellAnnoOrder, scale = F) +
  theme(panel.background = element_blank(), plot.background = element_blank())
## Loading required package: ggh4x
SpeDEG_Dotplot <- SpeDEG_Dotplot + scale_color_gradientn(colors = color.use, oob = scales::squish)
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
SpeDEG_Dotplot

Figure 5 E

WT_Mouse_SCPs <- subset(WT_Mouse, NewCelltype %in% c("(6) E Glia|SCP", "(11) PN Sustentacular"))

Genes = rev(c("Hey2", "Wnt6", "Col28a1", "Notch1","Ndrg1","Sfrp5",
              "Moxd1", "Sox5", "Gpr17", "Lmo4","Ngfr","Mpz",
              "Serpine2","Fabp7","Erbb3","Plp1","Sox2","Sox10"))
plot_list2 <- list()
for (g in Genes) {
  plot_list2[[g]] <- FeaturePlot(WT_Mouse_SCPs, features = g, order = T, alpha = 0.8, pt.size = 0.5, slot = "scale.data") + 
    theme(line = element_blank(), axis.text = element_blank(), axis.title = element_blank(), legend.position = "inside", legend.position.inside = c(0.1,0.7), plot.title = element_text(face = "bold.italic")) + 
    scale_color_gradientn(colors = c("dodgerblue", "grey90", "red"))
}
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
print(wrap_plots(plot_list2, ncol = 6) + theme(text = element_text(family = "Arial")))

Figure 5 F

OutDir <- "./DEG_GO/"
Whole_E_Glia_GO <- ComputeGO(GeneList = WT_Mouse_SpeDEGs$gene[which(WT_Mouse_SpeDEGs$cluster == "(6) E Glia|SCP")], ShowFigure = T, PlotFigure = T, SaveAsTable = T, FileName = "Whole_E_Glia_Enriched",
                             TopGO = 20, OrgDb = "Mouse",fromType = 'SYMBOL', toType = 'ENTREZID', Path = OutDir)
## Loading required package: clusterProfiler
## 
## Registered S3 method overwritten by 'ggtree':
##   method         from     
##   fortify.igraph ggnetwork
## clusterProfiler v4.10.1  For help: https://yulab-smu.top/biomedical-knowledge-mining-book/
## 
## If you use clusterProfiler in published research, please cite:
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation. 2021, 2(3):100141
## 
## Attaching package: 'clusterProfiler'
## The following object is masked from 'package:igraph':
## 
##     simplify
## The following object is masked from 'package:IRanges':
## 
##     slice
## The following object is masked from 'package:S4Vectors':
## 
##     rename
## The following object is masked from 'package:purrr':
## 
##     simplify
## The following object is masked from 'package:stats':
## 
##     filter
## Loading required package: org.Mm.eg.db
## Loading required package: AnnotationDbi
## 
## Attaching package: 'AnnotationDbi'
## The following object is masked from 'package:clusterProfiler':
## 
##     select
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## 'select()' returned 1:1 mapping between keys and columns
## Warning in bitr(GeneList, fromType = fromType, toType = toType, OrgDb =
## OrgDb_Used): 6.98% of input gene IDs are fail to map...
Whole_PN_Glia_GO <- ComputeGO(GeneList = WT_Mouse_SpeDEGs$gene[which(WT_Mouse_SpeDEGs$cluster == "(11) PN Sustentacular")], ShowFigure = T, PlotFigure = T, SaveAsTable = T, FileName = "Whole_PN_Glia_Enriched",
                              TopGO = 20, OrgDb = "Mouse",fromType = 'SYMBOL', toType = 'ENTREZID', Path = OutDir)
## 'select()' returned 1:1 mapping between keys and columns
## Warning in bitr(GeneList, fromType = fromType, toType = toType, OrgDb =
## OrgDb_Used): 2.22% of input gene IDs are fail to map...
E17_Glia_GO <- read.table(paste0(OutDir,"Whole_E_Glia_Enriched_GO_Sig.tsv"), header = T, sep = "\t")
PN_Glia_GO <- read.table(paste0(OutDir,"Whole_PN_Glia_Enriched_GO_Sig.tsv"), header = T, sep = "\t")

PN_Glia_GO_Spe <- PN_Glia_GO[-which(PN_Glia_GO$ID %in% E17_Glia_GO$ID),]

P_PNGlia_PathwayGO <- Bubble_Plot(pathway = PN_Glia_GO_Spe[grep("signaling", PN_Glia_GO_Spe$Description),][1:25,], "PN Glia Cell All Pathways")
P_PNGlia_PathwayGO + theme(text = element_text(family = "Arial"))

Figure 5 G

ucols_subset <- ucols[levels(CellChat_P90@idents)]
netVisual_individual(CellChat_P90, signaling = "NOTCH", pairLR.use = "DLK1_NOTCH1", layout = "circle", color.use = ucols_subset)

## [[1]]

Figure 5 H

ucols_subset <- ucols[levels(CellChat_AgedAdult@idents)]
netVisual_individual(CellChat_AgedAdult, signaling = "NOTCH", pairLR.use = "DLK1_NOTCH1", layout = "circle", color.use = ucols_subset)

## [[1]]

Figure 5 I

gg <- netVisual_bubble(object = cellChat_Merged, sources.use = c(1:3), targets.use = c(1:3),  comparison = c(1, 2, 3, 4, 5), angle.x = 90, signaling = c("WNT", "NOTCH"), return.data = T)
## Comparing communications on a merged object
SelectedCommunicationDF <- gg$communication

SelectedCommunicationDF_Modified <- SelectedCommunicationDF[grep("Wnt",SelectedCommunicationDF$interaction_name_2, invert = T),]
WNT_df <- SelectedCommunicationDF[grep("Wnt", SelectedCommunicationDF$interaction_name_2),]

WNT_df_Modified <- WNT_df %>% na.omit() %>% group_by(source.target, ligand) %>% reframe(
    prob = sum(prob, na.rm = TRUE),
    prob.original = sum(prob.original, na.rm = TRUE),
    pval = round(mean(pval, na.rm = TRUE)),
    interaction_name = paste(ligand, "- All Activated Receptors"),
    interaction_name_2 = paste(ligand, "- All Activated Receptors"),
    receptor = "All Activated Receptors",
    across(-c(prob, prob.original, pval, interaction_name, interaction_name_2), dplyr::first)) %>% ungroup()

WNT_df_Modified <- WNT_df_Modified[,colnames(SelectedCommunicationDF)]

SelectedCommunicationDF_Modified <- rbind(SelectedCommunicationDF_Modified, WNT_df_Modified)
g <- Figure5I_Quick(CellChatObj = cellChat_Merged, CommunicationDF = SelectedCommunicationDF_Modified)
g + theme(text = element_text(family = "Arial"))

Figure 5 J

Notch_Genes <- read.table("./Notch_Genes.csv", header = T, sep=",")
Notch_Genes$Gene <- DescTools::StrCap(Notch_Genes$Gene, method = "title")
## Registered S3 methods overwritten by 'proxy':
##   method               from    
##   print.registry_field registry
##   print.registry_entry registry
MyDotPlots(WT_Mouse_Medulla, Notch_Genes[which(Notch_Genes$Fine_Type != "Diff"),],GeneCol = 'Gene', Assay = "RNA", CelltypeCol = 'Fine_Type',Return = T,Draw = F, Group = "NewCelltype", CellTypeOrder = MedullaCellAnnoOrder, scale = F) +
  theme(panel.background = element_blank(), plot.background = element_blank()) + scale_color_gradientn(colors = color.use, oob = scales::squish)
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

Sup Figure 5 A - C

tmp1 <- DimPlot(WT_Mouse, group.by=c("samples"), reduction = 'umap', label = F) + theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), axis.line = element_blank()) + guides(color = guide_legend(ncol = 1))
tmp2 <- DimPlot(WT_Mouse, group.by=c("DetailedAge"), reduction = 'umap', label = F) + theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), axis.line = element_blank())
tmp3 <- DimPlot(WT_Mouse, group.by=c("seurat_clusters"), reduction = 'umap', label = T, repel = T) + theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), axis.line = element_blank()) + NoLegend()
tmp1 + tmp2 + tmp3

Sup Figure 5 D

MyMarkers <- read.table("./Mouse_AdrenalGland_Markers.txt", header = T, sep="\t")
MyMarkers$Markers <- DescTools::StrCap(MyMarkers$Markers, method = "title")
MyDotPlots(WT_Mouse, MyMarkers,GeneCol = 'Markers',Assay = "RNA", CelltypeCol = 'Celltype',Return = T,Draw = F, ClusterIdents = F, scale = F,
                      CellTypeOrder = c("6","11","10","8","0","2","1","5","18","7","14","12","9","3","4","13","16","15","19","17")) +
  theme(panel.background = element_blank(), plot.background = element_blank(), axis.text.x = element_text(vjust = 0.5, face = "italic"), axis.text.y = element_text(face = "plain"))

Sup Figure 5 E

ucols_subset <- ucols[levels(CellChat_AgedAdult@idents)]
netAnalysis_signalingRole_network(CellChat_AgedAdult, signaling = "NOTCH", width = 10, height = 3, font.size = 10, color.use = ucols_subset)

Sup Figure 5 F

ucols_subset <- ucols[levels(CellChat_AgedAdult@idents)]
netAnalysis_signalingRole_network(CellChat_AgedAdult, signaling = "WNT", width = 10, height = 3, font.size = 10, color.use = ucols_subset)

Sup Figure 5 G

ucols_subset <- ucols[levels(CellChat_AgedAdult@idents)]
netAnalysis_signalingRole_network(CellChat_AgedAdult, signaling = "BMP", width = 10, height = 3, font.size = 10, color.use = ucols_subset)

Sup Figure 5 H

layout(matrix(1:5, nrow = 1, byrow = TRUE))

ucols_subset <- ucols[levels(CellChat_E17@idents)]
p1 <- netVisual_individual(CellChat_E17, signaling = "NOTCH", pairLR.use = "DLK1_NOTCH1", layout = "circle", color.use = ucols_subset)

ucols_subset <- ucols[levels(CellChat_P5@idents)]
p2 <- netVisual_individual(CellChat_P5, signaling = "NOTCH", pairLR.use = "DLK1_NOTCH1", layout = "circle", color.use = ucols_subset)

ucols_subset <- ucols[levels(CellChat_P14@idents)]
p3 <- netVisual_individual(CellChat_P14, signaling = "NOTCH", pairLR.use = "DLK1_NOTCH1", layout = "circle", color.use = ucols_subset)

ucols_subset <- ucols[levels(CellChat_P90@idents)]
p4 <- netVisual_individual(CellChat_P90, signaling = "NOTCH", pairLR.use = "DLK1_NOTCH1", layout = "circle", color.use = ucols_subset)

ucols_subset <- ucols[levels(CellChat_AgedAdult@idents)]
p5 <- netVisual_individual(CellChat_AgedAdult, signaling = "NOTCH", pairLR.use = "DLK1_NOTCH1", layout = "circle", color.use = ucols_subset)

Sup Figure 5 I

layout(matrix(1:4, nrow = 1, byrow = TRUE))

ucols_subset <- ucols[levels(CellChat_P5_WNT@idents)]
netVisual_aggregate(CellChat_P5_WNT, signaling = "WNT", layout = "circle", color.use = ucols_subset) 
ucols_subset <- ucols[levels(CellChat_P14_WNT@idents)]
netVisual_aggregate(CellChat_P14_WNT, signaling = "WNT", layout = "circle", color.use = ucols_subset) 
ucols_subset <- ucols[levels(CellChat_P90_WNT@idents)]
netVisual_aggregate(CellChat_P90_WNT, signaling = "WNT", layout = "circle", color.use = ucols_subset) 
ucols_subset <- ucols[levels(CellChat_AgedAdult_WNT@idents)]
netVisual_aggregate(CellChat_AgedAdult_WNT, signaling = "WNT", layout = "circle", color.use = ucols_subset) 

Session Information

sessionInfo()
## R version 4.3.2 (2023-10-31)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: SUSE Linux Enterprise Server 15 SP5
## 
## Matrix products: default
## BLAS/LAPACK: /cfs/klemming/projects/supr/sllstore2017016/joey.d/Tools/anaconda3/envs/R432/lib/libopenblasp-r0.3.27.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Europe/Stockholm
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] org.Mm.eg.db_3.18.0         AnnotationDbi_1.64.1       
##  [3] clusterProfiler_4.10.1      ggh4x_0.2.8                
##  [5] showtext_0.9-7              showtextdb_3.0             
##  [7] sysfonts_0.8.9              CellChat_1.6.1             
##  [9] ggpubr_0.6.0                reshape2_1.4.4             
## [11] mclust_6.1.1                cluster_2.1.6              
## [13] cowplot_1.1.3               ggthemes_5.1.0             
## [15] ComplexHeatmap_2.18.0       gridGraphics_0.5-1         
## [17] gridExtra_2.3               extrafont_0.19             
## [19] ggrepel_0.9.5               patchwork_1.2.0            
## [21] igraph_2.0.3                network_1.19.0             
## [23] viridis_0.6.5               viridisLite_0.4.2          
## [25] SCopeLoomR_0.13.0           SingleCellExperiment_1.24.0
## [27] SummarizedExperiment_1.32.0 Biobase_2.62.0             
## [29] GenomicRanges_1.54.1        GenomeInfoDb_1.38.8        
## [31] IRanges_2.36.0              S4Vectors_0.40.2           
## [33] BiocGenerics_0.48.1         MatrixGenerics_1.14.0      
## [35] matrixStats_1.3.0           magrittr_2.0.3             
## [37] lubridate_1.9.3             forcats_1.0.0              
## [39] stringr_1.5.1               purrr_1.0.2                
## [41] readr_2.1.5                 tidyr_1.3.1                
## [43] tibble_3.2.1                ggplot2_3.5.1              
## [45] tidyverse_2.0.0             Seurat_5.0.0               
## [47] SeuratObject_5.0.2          sp_2.1-1                   
## [49] dplyr_1.1.4                 openxlsx_4.2.7             
## 
## loaded via a namespace (and not attached):
##   [1] gld_2.6.6               goftest_1.2-3           Biostrings_2.70.3      
##   [4] vctrs_0.6.5             spatstat.random_3.2-3   proxy_0.4-27           
##   [7] digest_0.6.35           png_0.1-8               shape_1.4.6.1          
##  [10] Exact_3.3               registry_0.5-1          deldir_2.0-4           
##  [13] parallelly_1.37.1       magick_2.8.3            MASS_7.3-60            
##  [16] httpuv_1.6.15           foreach_1.5.2           qvalue_2.34.0          
##  [19] withr_3.0.0             xfun_0.49               ggfun_0.1.5            
##  [22] survival_3.7-0          memoise_2.0.1           gson_0.1.0             
##  [25] systemfonts_1.1.0       ragg_1.3.2              tidytree_0.4.6         
##  [28] zoo_1.8-12              GlobalOptions_0.1.2     pbapply_1.7-2          
##  [31] KEGGREST_1.42.0         promises_1.3.0          httr_1.4.7             
##  [34] rstatix_0.7.2           globals_0.16.3          fitdistrplus_1.1-11    
##  [37] rstudioapi_0.16.0       miniUI_0.1.1.1          generics_0.1.3         
##  [40] DOSE_3.28.2             ggalluvial_0.12.5       zlibbioc_1.48.2        
##  [43] ggraph_2.2.1            polyclip_1.10-6         GenomeInfoDbData_1.2.11
##  [46] SparseArray_1.2.4       xtable_1.8-4            doParallel_1.0.17      
##  [49] evaluate_0.24.0         S4Arrays_1.2.1          hms_1.1.3              
##  [52] bookdown_0.41           irlba_2.3.5.1           colorspace_2.1-0       
##  [55] hdf5r_1.3.11            ggnetwork_0.5.13        ROCR_1.0-11            
##  [58] readxl_1.4.3            reticulate_1.38.0       spatstat.data_3.0-4    
##  [61] lmtest_0.9-40           ggtree_3.10.1           later_1.3.2            
##  [64] lattice_0.22-6          spatstat.geom_3.2-9     NMF_0.28               
##  [67] future.apply_1.11.2     shadowtext_0.1.4        scattermore_1.2        
##  [70] RcppAnnoy_0.0.22        class_7.3-22            pillar_1.9.0           
##  [73] nlme_3.1-165            iterators_1.0.14        sna_2.8                
##  [76] gridBase_0.4-7          compiler_4.3.2          RSpectra_0.16-1        
##  [79] stringi_1.8.4           DescTools_0.99.58       tensor_1.5             
##  [82] plyr_1.8.9              crayon_1.5.2            abind_1.4-5            
##  [85] haven_2.5.4             graphlayouts_1.2.0      bit_4.0.5              
##  [88] rootSolve_1.8.2.4       fastmatch_1.1-4         textshaping_0.3.7      
##  [91] codetools_0.2-20        bslib_0.7.0             e1071_1.7-14           
##  [94] lmom_3.2                GetoptLong_1.0.5        plotly_4.10.4          
##  [97] mime_0.12               splines_4.3.2           circlize_0.4.16        
## [100] Rcpp_1.0.12             fastDummies_1.7.3       HDO.db_0.99.1          
## [103] cellranger_1.1.0        Rttf2pt1_1.3.12         knitr_1.47             
## [106] blob_1.2.4              utf8_1.2.4              clue_0.3-65            
## [109] fs_1.6.4                listenv_0.9.1           expm_1.0-0             
## [112] ggsignif_0.6.4          ggplotify_0.1.2         Matrix_1.6-5           
## [115] statmod_1.5.0           tzdb_0.4.0              svglite_2.1.3          
## [118] tweenr_2.0.3            pkgconfig_2.0.3         tools_4.3.2            
## [121] cachem_1.1.0            RSQLite_2.3.7           DBI_1.2.3              
## [124] fastmap_1.2.0           rmarkdown_2.27          scales_1.3.0           
## [127] ica_1.0-3               broom_1.0.6             sass_0.4.9             
## [130] coda_0.19-4.1           FNN_1.1.4               BiocManager_1.30.23    
## [133] dotCall64_1.1-1         carData_3.0-5           RANN_2.6.1             
## [136] farver_2.1.2            scatterpie_0.2.4        tidygraph_1.3.1        
## [139] yaml_2.3.8              cli_3.6.2               leiden_0.4.3.1         
## [142] lifecycle_1.0.4         uwot_0.2.2              mvtnorm_1.2-5          
## [145] presto_1.0.0            backports_1.5.0         BiocParallel_1.36.0    
## [148] timechange_0.3.0        gtable_0.3.5            rjson_0.2.21           
## [151] ggridges_0.5.6          progressr_0.14.0        ape_5.8                
## [154] parallel_4.3.2          limma_3.58.1            jsonlite_1.8.8         
## [157] RcppHNSW_0.6.0          bitops_1.0-7            bit64_4.0.5            
## [160] Rtsne_0.17              yulab.utils_0.1.8       spatstat.utils_3.0-5   
## [163] BiocNeighbors_1.20.2    zip_2.3.1               jquerylib_0.1.4        
## [166] highr_0.11              GOSemSim_2.28.1         lazyeval_0.2.2         
## [169] shiny_1.8.1.1           htmltools_0.5.8.1       enrichplot_1.22.0      
## [172] GO.db_3.18.0            sctransform_0.4.1       glue_1.7.0             
## [175] spam_2.10-0             XVector_0.42.0          RCurl_1.98-1.14        
## [178] treeio_1.26.0           boot_1.3-30             extrafontdb_1.0        
## [181] R6_2.5.1                labeling_0.4.3          rngtools_1.5.2         
## [184] aplot_0.2.3             statnet.common_4.10.0   DelayedArray_0.28.0    
## [187] tidyselect_1.2.1        ggforce_0.4.2           car_3.1-2              
## [190] future_1.33.2           munsell_0.5.1           KernSmooth_2.23-24     
## [193] data.table_1.15.4       htmlwidgets_1.6.4       fgsea_1.28.0           
## [196] RColorBrewer_1.1-3      rlang_1.1.4             spatstat.sparse_3.0-3  
## [199] spatstat.explore_3.2-7  Cairo_1.6-2             fansi_1.0.6